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O ' Based on a coherent state representation of noise operator and an ensemble averaging procedure 

we have recently developed [Phys. Rev. E 65, 021109 (2002); ibid. 051106 (2002)] a scheme for 
quantum Brownian motion to derive the equations for time evolution of true probability distribution 
I ' functions in c-number phase space. We extend the treatment to develop a numerical method for 

generation of c-number noise with arbitrary correlation and strength at any temperature, along with 
■j-J , the solution of the associated generalized quantum Langevin equation. The method is illustrated 

with the help of a calculation of quantum mean first passage time in a cubic potential to demonstrate 
' quantum Kramers turnover and quantum Arrhenius plot. 



PACS numbers: 05.40-a, 82.20.- 



I. INTRODUCTION 



A system coupled to its environment is the standard paradigm for quantum theory of Brownian 
motioniiSi&iSiSiLSiiiiifliiLiS. Its overwhelming success in the treatment of various phenomena in quantum optics, 
transport processes in Josephson junction, coherence effects and macroscopic quantum tunneling in condensed matter 
0^ ' physics, electron transfer in large molecules, thermal activation processes in chemical reactions is now well known 
and forms a large body of current literature. While the early development of quantum optics initiated in the sixties 
and seventies was based on density operator, semigroup, noise operator or master equation methods primarily within 
weak-coupling and Markov approximations, path integral approach to quantum Brownian motion attracted wide at- 
tention in the early eighties. Although this development had widened the scope of condensed matter and chemical 
physics significantly, so far as the large coupling between the system and the heat bath and large correlation times of 
the noise processes are concerned several problems still need to be addressed. First, a search for quantum analogue 
of Kramers' equation for a nonlinear system which describes quantum Brownian motion in phase space had remained 
elusive and at best resulted in equation of motion which contains higher (than second) derivatives of probability dis- 
tribution functions 2,13 whose positive definiteness is never guaranteed. As a result these quasi-probability distribution 
functions often become singular or negative in the full quantum domain. Second, although long correlation times and 
large coupling constants are treated non-perturbatively formally in an exact manner by path integrals, their analytic 
evaluation is generally completed within semi-classical schemesi^iii. This results in situations where the theory fails 
to retain its validity in the vacuum limit. Third, although the numerical techniques based on the path integral Monte 
Carlo are very successful in the treatment of equilibrium properties it is often difficult to implement in a dynamical 
scheme because of the well-known sign problem associated with varying phases in a quantum evolution over paths due 
to the oscillating nature of the real time propagatoriiiiS*^. Fourth, the treatment of non-Markovian dynamics within 
a quantum formulation as retardation effect in a memory functional is very difficult, if not impossible, to handle even 
numericallyi°*i£. 

Keeping in view of these difficulties it is worthwhile to ask how to extend classical theory of Brownian motion 
to quantum domain for arbitrary friction and temperature. Based on a coherent state representation of the noise 
operator and an ensemble averaging procedure we have recentl}>i^*i^S4ii developed a scheme for quantum Brownian 
motion in terms of a c-number generalized Langevin equation. The approach allows us to use classical methods of non- 
Markovian dynamics to derive exact generalized quantum Kramers' equation and other relevant quantum analogues 
of classical diffusion, Fokker-Planck and Smoluchowski equations. The object of the present analysis is to extend 
the treatment (?) to develop a numerical method for generation of quantum noise as classical c-numbers («) and to 
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solve the associated generalized quantum Langevin equation. The numerical method valid for arbitrary temperature, 
damping strength and noise correlation, is exemplified by an application to calculation of quantum mean first passage 
time in a cubic potential. The approach is classical-looking in form, and is independent of path integral quantum 
Monte Carlo techniques and takes care of quantum effect order by order with simplicity and accuracy to a high degree. 

The outlay of the paper is as follows: We start with a brief review of c- number representation of quantum noise and 
the associated quantum Langevin equation, as formulated recently by us, in Sec. II. This is followed by a numerical 
method for generation of quantum noise and solution of the Langevin equation in Sec. III. In Sec. IV we illustrate the 
procedure by calculating quantum mean first passage time in a cubic potential with an emphasis on turn-over features 
and Arrhenius plot in thermal activated processes assisted by tunneling. 



II. A COHERENT STATE REPRESENTATION OF QUANTUM NOISE AND QUANTUM LANGEVIN 

EQUATION IN C-NUMBERS 

We consider a particle of unit mass coupled to a medium comprised of a set of harmonic oscillators with frequency 
u>i. This is described by the following Hamiltonian: 
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Here x and p are co-ordinate and momentum operators of the particle and the set {<lj,Pj} is the set of co-ordinate and 
momentum operators for the reservoir oscillators coupled linearly to the system through their coupling coefficients 
Kj. The potential V(x) is due to the external force field for the Brownian particle. The co-ordinate and momentum 
operators follow the usual commutation relations [x,p] = ih and [qj,pj] — ihSij. Note that in writing down the 
Hamiltonian no rotating wave approximation has been used. 

Eliminating the reservoir degrees of freedom in the usual way^Li^ we obtain the operator Langevin equation for 
the particle, 



x(t) + / dt'~/(t - t')x(t') + V'(x) = F(t) , (2) 
Jo 

where the noise operator F(t) and the memory kernel j(t) are given by 

F(t) = [{<7j(°) - *(°)} K i cosujjt+Pj^K 1 / 2 smujjt (3) 

3 

and 

7(*) = J^KjCOSujjt , (4) 

3 

with Kj — u)j ( masses have been assumed to be unity ). 

The Eq.J2l is an exact quantized operator Langevin equation^ for which the noise properties of F(t) can be defined 
using a suitable initial canonical distribution of the bath co-ordinates and momenta. Our aim here is to replace 
it by an equivalent QGLE in c-numbers. To achieve this and to address the problem of quantum non-Markovian 
dynamics in terms of a true probabilistic description we proceed in two steps^i^iSiii. We first carry out the quantum 
mechanical average of Eq.J2J) 

(£(*)) + f dt' 7 (t-t'){x(t')) + (V'(x)) = (F(t)) (5) 



n 



where the average (...) is taken over the initial product separable quantum states of the particle and the bath oscillators 
at t = 0, |</>){|ai)|a2) ■ ■ • |ajv)}- Here \4>) denotes any arbitrary initial state of the particle and \ati) corresponds to 
the initial coherent state of the z-th bath oscillator. is given by |a<i) = exp(— |ai| 2 /2) Y^ i= o( a T I *J™J-)\ n i) > a i 
being expressed in terms of the mean values of the co-ordinate and momentum of the i-th oscillator, (qi(0)) — (x(0)) = 
{\fhj2iOi)(ai + a*) and (pi(0)) — i\ffojjjj2(a* — a*), respectively. It is important to note that (F(t)) of Eq.© is a 
classical-like noise term which, in general, is a non-zero number because of the quantum mechanical averaging over the 
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co-ordinate and momentum operators of the bath oscillators with respect to the initial coherent states and arbitrary 
initial state of the particle and is given by 

(F(t)) = [{<&(0)> - (4(0))} Kj cos^t + (fr (0))«} /a sinwjt] . (6) 
j 

It is convenient to rewrite the c-number equation as follows; 

(x(t)) + f df 7 (* - 0(*(0) + (*"(*)) = F(t) (7) 



where we let the quantum mechanical mean value (F(t)} = f(t). We now turn to the second averaging. To realize F(t) 
as an effective c-number noise we now assume that the momenta (pj(0)) and the shifted co-ordinates {(qj(0)) — (i(0))} 
of the bath oscillators are distributed according to a canonical distribution of Gaussian forms as 

^■ = AAexp{-[fe^ 2 + ^;^W)"(x(0))} 2 ]| 
[ («j + 2) J 

so that for any quantum mechanical mean value Oj((pj(0)) , ((<?j(0)) — (x(0))) the statistical average (. . .)g is 



(Oi)s = J Oj((Pj(o)),{(?,-(o))-<i(o))}) 

x^(fe(0)),{fe(0))-(x(0))}) 

xd(fc(0))d{<&(0)) - (x(0))} . (9) 

Here indicates the average thermal photon number of the j-th oscillator at temperature T and fij — 
l/[exp (hujj /ksT) — 1] and Af is the normalization constant. 

The distribution © and the definition of statistical average (jHJl imply that F(t) must satisfy 

(F(t)) s = (10) 

and 

(F(t)F(t')) s = i Yl ( coth ^) cos "i(t - t') (11) 

That is c-number noise F(t) is such that it is zero-centered and satisfies the standard fluctuation-dissipation relation 
(FDR) as known in the literature in terms of the quantum statistical average of the noise operators. We emphasize 
that the method of ensemble averaging over c-number noise F(t) is distinctly different from traditional methods of 
quantum statistical averaging and also from the method of Glauber-Sudarshan P-distribution. We discuss this in 
Appendix A. 

We now add the force term V'((x)) on both sides of Eq.Q and rearrange it to obtain 

X = P (12a) 



P = - dt' 1 {t-t')P{t')-V'{X) + F{t) + Q{t) (12b) 
Jo 

where we put (x(t)) = X(t) for notational convenience and 



Q(t) = V'((x)) - (V'(x)) (13) 

represents the quantum correction due to system degrees of freedom. Since Q(t) is a fluctuation term Ea. (|12b)l offers 
a simple interpretation. This implies that the QGLE is governed by a c-number quantum noise F(t) originating 
from the heat bath characterized by the properties (|1(J|) and (|ll|l and a quantum fluctuation term characteristic of 
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the non-linearity of the potential. To proceed further we need a recipe for the calculation of Q(t). This has been 
discussed earlier in several contextsi^*iSiiSiiLi2iiS. For the present analysis we summarize it as follows. 
Referring to the quantum nature of the system in the Heisenberg picture, one may write. 

x(t) = (x) + Sx (14a) 
p(t) = (p)+Sp (14b) 

where (x) and (p) are the quantum- mechanical averages and Sx, dp are the operators. By construction (Sx) and (dp) 
are zero and [Sx, Sp] = ih. Using Eas. (|14a|l and 114b|) in (V'(x)) and a Taylor series expansion around (x) it is possible 
to express Q(t) as 

Q(*) = -E^ F(n+1) w^ n (*)) ( 15 ) 

n>2 

Here V n (X) is the n-th derivative of the potential V(X). To the second order Q(t) is given by Q(t) = -\V" '(X)(8x 2 ) 

where X(t) and (Sx 2 (t)) can be obtained as explicit functions of time by solving the following set of coupled equations 
(HSJ-CB) along with Eos. ((Hal) and fTSEfr : 

d(Sx 2 )/dt = (SxSp + SpSx) (16) 

d(SxSp + SpSx)/dt = 2(Sp 2 ) -2V"(X)(Sx 2 ) (17) 

d(Sf)/dt = -V"{X)(SxSp + SpSx) (18) 

We now make some pertinent remarks on the above equations for quantum corrections and the related approxima- 
tions and their successive improvements. The set of equations <|16fl - <|18[) f° r calculation of quantum correction term 
to the lowest order on based on 8 , the equations Sx = Sp and Sp = —V"(X)Sx, when the system-bath correlation is 
taken into account through X in V"(X) which is obtained as a solution to Eas. (|12a|l and (|12b|) . The system-reservoir 
correlation is included at the non-Markovian level in the dynamics involving quantum correction terms. This implies 
that we have neglected the two time correlation of the system like (Sp(i)Sp(t 1 )) . The next improvement would be to 
include the contribution due to these terms at least in the Markovian level to the quantum correction equations so 
that the right hand sides of Eas. i|17fl and i|18|) may be modified to add — Y(Sx5p + SpSx) and — 2Y(Sp 2 ), respectively, 
r is the Markovian limit of dissipation. The quantum correction dynamics is thus taken into consideration at the 
Markovian level through Y term and at the non-Markovian level through V"(X) term while the mean values follow 
non-Markovian dynamics. Being exponential in nature the contribution of the former terms are felt only at long times 
when the quantum fluctuations in the low order has less significant role to play with. Furthermore, it is now well 
establishedSiMiSiSi that curvature of the potential V"(X) controls the initial growth of the quantum fluctuations. 
Since X(t) is stochastic V"(X) imparts stochasticity in Sx = —V"(X)Sx, which makes the average dynamics of (Sx 2 ) 
exponentially diverging 38 . We thus emphasize that the V"(X) term dominantly controls the evolution of quantum 
fluctuations to a good degree of accuracy. In the overdamped limit, however, it is easy to incorporate quantum cor- 
rection to very high orders. To this end one may start from equations Sx = Sp and Sp — — YSp— V"(X)Sx. Neglecting 
the inertial term (SX) in this limit we may write —YSx = V"(X)5x to yield the solutions (Sx n (t)) for arbitrary n in 
the expression for Q(t). This has been discussed in greater detail inii in the of quantum Smoluchowskii equation. For 
a more accurate estimate of quantum corrections one may use higher order terms in Q(t) which can be calculated by 
solving higher order (say, twelve equations for fourth order) coupled equations as derived in Appendix B. Under very 
special circumstance, it has been possible to include quantum effects to all orders^. 

The above method of c-number representation of quantum Langevin equation has been utilized recently by us to 
derive exact quantum Kramers' equationi£*i£ in terms of the probability distribution function for calculation of escape 
rate from a metastable well and to address several other issues. In what follows we present a method for generation 
of c-number noise for numerical simulation of quantum Langevin equation i|12a|) and l|12b|l in the next section. 

III. A NUMERICAL METHOD FOR GENERATION OF QUANTUM NOISE 

Eq. i|ll|) . the fluctuation-dissipation relation is the key element for generation of quantum noise. The left hand side 
of this equation suggests that F(t) is a 'classical' random number and (F(t)F(t')) s is a correlation function which 
is classical in form but quantum-mechanical in its content. That is F(t) is a c-number Langevin force such that 
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Ea. (|12bjl (along with H12afl ) represents a stochastic differential equation. Since one of the most commonly occuring 
noise processes is exponentially correlated colour noise, we develop the scheme first for this process. This is based 
on a simple idea that c-number quantum noise may be viewed as a superposition of several Ornstein-Ulhenbeck noise 
processes. It may be noted that Imamoglu^i, a few years ago, has produced a scheme for approximating a class of 
reservoir spectral functions by a finite superposition of Lorentzian functions with positive coefficients corresponding 
to a set of fictitous harmonic oscillator modes to incorporate non-Markovian effects into the dynamics. The aim of the 
work21 is to develop a non-Markovian dynamics in terms of a master equation corresponding to the system-reservoir 
interaction for the expanded model after the tracing operation. The method is an extension of the stochastic wave 
function method in quantum optics. The present superposition method is reminiscent of this approach in the time 
domain. We however emphasize two important points. Because of rotating wave and Born approximations the master 
equation is valid for weak system- reservoir coupling i.e., weak damping regime. The present method on the other hand 
uses a superposition of noises in the time domain and the quantum Langevin equation is valid for arbitrary coupling 
since no weak coupling approximation or rotating wave approximation has been used in the treatment. Secondly, 
while the stochastic wave function approach as well as the density operator equation have no classical analogue, the 
present scheme makes use of c-numbers to implement the classical numerical method quite freely to solve quantum 
Langevin equation. 

A. SUPERPOSITION METHOD 

We proceed in several steps as follows; 

1. Step I 

We first denote the correlation function c(t — t') as 

(F(t)F(t')) s = c(t-t') (19) 
and in the continuum limit c(t — t') is given by 

c(t-t') = ~ diJK(u)p(u)hw coth(— — - ) cos w(i-t') (20) 
2 J 2k B T 

where we have introduced the density function p(uj) for the bath modes. In principle an a priori knowledge of k(lo)p(lo) 
is essential and primary factor for determining the evolution of stochastic dynamics governed by Eas. (|12a|) and l|12b(l . 
We assume a Lorentzian distribution of modes so that 

k ^ = It^w (21) 

where T and r c are the dissipation in the Markovian limit and correlation time, respectively. These two quantities 
and the temperature T are the three input parameters for our problem. Ea. (|21|l when used in |0J in the continuum 
limit yields an exponential memory kernel j(t) = (Y/t c ) exp(— |t|/r c ). For a given set of these parameters we first 
numerically evaluate the integral (|20|l as a function of the time (t — t'). The typical profile of such a correlation is 
shown in Fig. 1(a) for three different temperatures (solid line) and for fixed T = 1 and r c = 1. Similarly we present the 
numerically calculated correlation functions for different values of Y for a fixed temperature T = 0.0 and correlation 
time t c = 1.0; (Fig 1(b), solid line) and for different values of r c (for fixed temperature T = 0.0 and dissipation 
r = 1.0; (Fig 1(c), solid line). 

2. Step II 

In the next step we numerically fit the correlation function with a superposition of several exponentials as 



(22) 
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The set Di and r,, the constant parameters are thus known. In Fig.l (a-c) we exhibit these fitted curves (dotted 
lines) against the corresponding solid curves calculated by numerical integration of Ea. (|2U[) . It may be checked that 
at high temperatures the required fitting curve is in general a single exponential, whereas at very low temperatures 
down to absolute zero biexponential fitting is an excellent description. At the intermediate temperatures one has to 
use even triexponential fitting, in general. 

3. Step-III 

The numerical agreement between the two sets of curves based on 112011 and ill'L't in each of the Figs. 1 (a-c) suggests 
that one may generate a set of exponentially correlated colour noise variables rji according to 

m = -- + -m (23) 

Ti Ti 

where 



&(*)} = (24) 
(6(0)6(t)> = 2AM(t) (» = 1,2,3...) (25) 

in which £j(i) is a Gaussian white noise obeying (|24|) and l)25[l: t% and Di are being determined in Step-II from 
numerical fit. The quantum noise rji is thus an Ornstein-Ulhenbeck process with properties 

(Vi(t)) = (26) 
{ m {t) m {t')) = %^e X p(±^T) (27) 

Clearly Ti and Di are the correlation time and the strength of colour noise variables rji, respectively. It is also 
important to point out that they have nothing to do with true correlation time t c and strength T of the noise F(t) as 
expressed in (j?T)l and (|2"U|l . 



4. Step-IV 

The c-number quantum noise F(t) due to the heat bath is therefore given by 



F(t)= Vi (28) 

i=l,2.. 

Eq.l|28p implies that F(t) can be realized as a sum of several Ornstein-Ulhenbeck noises rji satisfying (F(i)F(t'))s — 

Having obtained the scheme for generation of c-number quantum noise F(t) we now proceed to solve the quantum 
Langevin equations H12al) and l|12bl) along with Ecis. Hl(jfl - i|18fl . To this end the first step is the integration22i2i of the 
stochastic differential equation J22J). We note that to generate Gaussian white noise £i(t) of from two random 
numbers which are uniformly distributed on the unit interval, Box-Mueller algorithm may be used. To be specific we 
start the simulation with integration of (|23|l to generate an initial value of rji as follows; 



m 
n 

Vi 



— random number 
= random number 



2D, 



■ ln(m) 



1/2 



cos(27rn) 



(29) 
(30) 

(31) 



We then set Ei = exp(— At/rj), where At is the integration step length for (|23fl . Following the integral algorithm of 
Fox et a@&, we now generate the next step rji\t+At from rji according to the following lines; 
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a = random number 
b = random number 

~(l-£?)ln(a) 



1/2 



cos(27r&) 



(32) 
(33) 

(34) 



and 



ViU+At = ViEi + hi (35) 

After Ij35(l the algorithm is repeated back to (|32J) and continued for a sufficient time evolution. 

In the next step we couple the above-mentioned integral algorithm for solving (|23|l to standard fourth order Runge- 
Kutta (or other) algorithm for the solution of Eas. i|12a|l . <|12b[) and H16fl - I|18|) . For an accurate and fast algorithm 
r\- integration (|35(l may be carried out at every At/2 step, while for Runge-Kutta integration of Eas. p2a(l . I|12b(l and 
(|16[) - <|18[1 step size At is a good choice. With At in the range ~ 0.001 — 0.01 the overall method gives excellent results 
for the solution of quantum Langevin equation. 



B. SPECTRAL METHOD 



Our method of superposition of Ornstein-Ulhenbeck noises for generation of quantum noise as discussed above pri- 
marily pertains to situations where the correlation function can be described by a superposition of several exponentials. 
Although in overwhelming cases of applications in condensed matter and chemical physics at low-temperature where 
in general non-Ohmic conditions prevail we encounter such situations, there are cases where the noise correlation 
exhibits power law behaviour or is Gaussian in nature. To extend the present scheme we now show how the spec- 
tral method can be implemented for generation of quantum noise. The method had been used earlier by several 
group922iSi2ii2^ in the context of classical long and short range temporally correlated noises. We give here a brief 
outline for its extension to quantum noise. 

The starting point of our analysis is the knowledge of quantum correlation function c(t — t 1 ) as given in J5DJ. 
Again this depends on k(uj) p(to) , which (for example, for short range Gaussian correlated noise may be taken as 
2rexp(— aj 2 r,?/4), T and r c being the strength and the correlation time, respectively) must be known a priori. Having 
known c(t — t') as a function of time we first carry out its Fourier transform 

(F(t)F(t')) s = c(t-t') 

C{uj) = J exp(~iujt)c{t)dt (36) 

by discretizing the time in N — 2" intervals of size At (such that At is shorter than any other characteristic time of 
the system). The result is a string of discrete numbers C{uo n ). In discrete Fourier space the noise correlation is given 
by 

(F(Lu n )F(Lo' m ))s = NAtC(uj n )6 m+ni0 (37) 
Therefore the noise in the Fourier domain can be constructed simply as 



F(w«) = y/NAtC(LJ n )a(u n ) (38) 
where i = 0, 1, 2, N; F(lj ) = F{w n ); w„ = j^- 

Here a n are the Gaussian complex random numbers with zero mean and anti-delta correlated noise, such that a n = 
and b n (the real and imaginary parts of a n ) are Gaussian random numbers with zero mean value and 

variance 1/2. 
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(ch) = (b n )=0; (39) 
K) = K) = \] where n f (40) 
and a\ = 1 (41) 

Box-Mueller algorithm may again be used to generate the set of random numbers a n , b n (and consequently a n ) 
according to the lines; 

p n = random number (42) 

q n = random number (43) 

a n = (-lnp n ) 1/2 cos27T(7„ (44) 

and similarly for b n . The discrete Fourier transform of the string F(u} n ) is then numerically calculated by the Fast 
Fourier Transform algorithm to generate the quantum random numbers F(t n ). The recipe can be implemented for 
solving Eas. H12al) and (|12b(l and l|16l) - (|18[) as discussed in the earlier subsection UFA. 



IV. APPLICATION: QUANTUM MEAN FIRST PASSAGE TIME CALCULATIONS 

We now present the numerical results on c-number quantum Langevin equation for a cubic potential. The mean first 
passage time whose inverse represents the quantum rate of crossing the barrier is computed for different temperatures, 
correlation times and damping constants. The notion of mean in mean first passage time calculation needs some 
clarification. We are concerned here with evolution of quantum mean values q{= (q)) and p{= (p)) which can be 
simultaneously determined unlike their operator counterparts q and p, respectively. The quantum fluctuation around 
the mean path are however taken into account by following the equations of quantum corrections. Since the system is 
coupled to heat bath it is also necessary to take statistical 'mean' over several thousands of realizations of trajectories 
of the quantum mean values. One of the major aims of this numerical simulation is to show that the method described 
is subsection UFA can be implemented very easily to describe both thermally activated processes and tunneling within 
an unified scheme. 

To begin with we consider a potential of the form 



V(x) — ——Ax 3 + Bx 2 (45) 

where A and B are the constants. The other input parameters for our calculations fEas. H20[) and $21$ ) are temperature 
T, strength of noise correlation T and correlation time r c . For the potential V(x), Q(t) is given by a single quantum 
correction term Q = A(Sx 2 ) according to Ea. (|15|l . Theoretical consideration demands that the activation barrier AE 
must be much larger than the thermal noise, i.e, AE 3> fcsT, where AE denotes the height of the potential barrier. 
In order to ensure the stability of the algorithm we have kept At/r c -C 1 where At is the integration step size. For 
a given parameter set wc have computed the quantum first passage time for the quantum Brownian particle in the 
cubic potential starting from the bottom of the metastable well at x = to reach the maximum at x = 2B /A for a 
single realization of a stochastic path. The quantum mean first passage time (QMFPT) is calculated by averaging 
over the first passage times over 5000 trajectories. We have used h = ks = 1 for the entire treatment. 

We now present the results for the quantum mean first passage times obtained by the simulation of the quantum 
damped cubic potential (|45|l driven by quantum colour noise [Eas. (|12a|) . l|12bjl . (|16fl - l|18|) and (|23|l ]. Our results are 
shown for the rate 1/QMFPT as a function of T for three different temperatures T = 3, 1 and 0.0 in Fig.2(a-c). We 
have kept the parameter AE — 12 with A = 0.5 and r c very small. It is observed that at T — 3 with increasing T, 
1/QMFPT undergoes a turnover from an inverse trend (1/MFPT cx 1/T) to a linear behaviour (1/QMFPT cx T) in 
accordance with Kramers' classic result subsequently supported by many thers2£*22i22422i2L2S. The 1/QMFPT has 
a maximum in the intermediate range. As the temperature is lowered to T = 1 and T = the turnover persists with 
a little change of shape near the maxima. The Fig. 2c represents the typical vacuum field assisted tunneling since the 
thermal activated process is almost absent at T ~ 0. 

In Fig. 3 we plot the rate 1/QMFPT vs. r c for two different temperatures, with A = 1.0 and AE = 4.0 and T = 1.0. 
It is observed that as the temperature is lowered to T = the decrease of 1/QMFPT exhibits two distinct time scales. 
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In Fig.4(a-b) we exhibit Arrhenius plot, i.e, QMFPT -1 vs. 1/T for several values of T. We have fixed AE = 10.0 
with A = 1.5 and r c = 1.0. The physical situation corresponding to non-Ohmic underdamped to overdamped 
regions has been explored. It is apparent that when the temperature is high InQMFPT -1 vs. 1/T is typically linear 
corresponding to the century old Arrhenius classical result. However, when the temperature is lowered the quantum 
effects and T 2 dependence is observed. This result is in complete conformity with those reported earlier in connection 
with the theoretical as well as experimental studies in mixed oxide systems and gas kinetic o 1 ! 1 '^ . The variation of 
Arrhenius profiles with T is also indicative of the turnover phenomenon as discussed earlier. 



V. CONCLUSION 



In this paper we have proposed a simple numerical method for generation of quantum noise and for solution of 
the associated generalized quantum c-number Langevin equation. The method is based primarily on two aspects. 
First, it has been shown that it is possible to realize quantum noise (in a quantum Langevin dynamics) which 
satisfies quantum fluctuation-dissipation relation as c-numbers rather than as operators whose ordering and quantum 
statistical averaging procedure is markedly different from those followed in the present scheme. Second, the quantum 
correction terms which appear due to the non-linearity of the potential of the system, can be generated systematically 
order by order, so that the coupled set of equations for the stochastic dynamics of the system and the correction 
terms can be solved easily with sufficient degree of accuracy. The method of simulation is equipped to deal with 
arbitrary correlation time, temperature and damping strength. The mapping of quantum stochastic dynamics as a 
problem of classical simulation as presented here allows us to implement classical methods in quantum simulations. 
The method is much simpler to handle compared to path integral quantum Monte-Carlo as employed earlier in the 
problem of barrier crossing dynamics^. We have shown that our simulation results on quantum mean first passage 
time for a particle in a cubic metastable potential agree well with the results obtained earlier. We hope that with the 
prescription given, many problems where the time correlation does not follow a typical prescribed dynamics can be 
simplified within the framework of this simulation. 
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APPENDIX A: FLUCTUATION-DISSIPATION RELATION 



In this section we derive quantum fluctuation-dissipation relation in three different ways. 



Method-I: Quantum-Statistical Averaging 

To make the discussion self-contained we first establish that the quantum-statistical average value of a two-time 
correlation function of the random force operator F(t) using standard usage of quantum-statistical averaging (...)qs 
can be written as 



^(F(t)F(f)+F(t')F(t)) QS =^^H (coth^) co8 Wj -(t-f) 



(Al) 



where F(t) is defined by Eq.©. Starting from a canonical distribution of bath oscillators at t = with respect to 
the free oscillator Hamiltonian Hq, with a shift of origin q'j(0) = qj(0) — x(0), where 

^ = Ef^ + ^f(0)), (A2) 
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we construct the expectation value of an arbitrary operator O as 



(6) 



QS 



Tr C> e -go/^ T 
Tr e-Ha/k B T 



(A3) 



As usual Trace implies both quantum mechanical averaging with respect to number states of the oscillators 
n^_ {|ni)}n^_ {|n2)}...n^°_ {|njv)} and the system state \<f>) followed by Boltzmann statistical averaging. Making 
use of the definition of the quantum statistical average Eq. I|A3|) we first derive 



(OTft(o)) 



QS 



(Pi(O)Pfc(O)) 



QS 



Ojk- 



hcot\\{uj,jh/2k B T) 



2ujj 



(A4) 



and 



(q'(0)p k (0))QS = -{pk(0)q'(0)) Q s = Mk 



(A5) 



Since the expectation value of an odd number of factors of <7j(0) and Pj(0) vanishes and the expectation value of 
an even number of factors is the sum of products of the pairwise expectation values with the order of the factors 
preserved we note that Gaussian property of the distribution of <?j(0) and Pj(0) is implied. 

The relations (|A4|) and (| A5|) directly yield (|A1|) . On the other hand quantum mechanical averaging with respect 
to the number states results in 



(F(t)) QS = 



(A6) 



Method-II: Glauber-Sudarshan P-function 

In this method we first calculate the quantum mechanical average of l/2[F(t)F(t') + F(t')F(t)} with respect to 
the coherent states of the bath oscillators |ai)|a2)...|ajv) and system state \(j>) at t = where \ai) has been defined 
in the text. \ai) is the eigenstate of the annihilation operator a, such that a,i\ai) = ai\a.i) where <^(0) and pi(0) 
are expressed as usual as <7i(0) = ^/h/2uji[a\{Q) + a,i(0)} and pi(0) — i\J hu)i/2[a\{0) — fij(Q)]. We thus note that the 
coherent state averaging (...) leads to (we have discarded the cross terms like aja% since they will vanish after the 
statistical averaging that follows in the next step) , 



hF(t)F(t') + F(t>)F(t)) = ]T 



Kj- — - [a* 2 + 2a* a j + 1 + ayj cos ujjt cosujjt' 



+ —KjLUj (2a* 2 — 2a 2 ^j sinuij(t — t') 



2a* a j — 1 + a ,■ J sinw^i sinujjt' 



(A7) 



We now introduce Glauber-Sudarshan distribution known as P-function for the statistical averaging of the quantum 
mechanical expectation values, 



P GS (a,(0),a*(0)) = A^exp 



A^exp 



fe(0)) 2 +t^(0))^ 
2hio jfij 



(A8) 



where N is the normalization constant. Noting that with the distribution (|A8|) the normalized statistical averages 
(— )gs are given by, 



{dj)GS = (a* )gs = and (aja* k ) GS = nj5 jk 



(A9) 
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With relations i|A9|l , the statistical average (with respect to the P- function) of the quantum mechanical mean value 
(|A7(I leads us to 

\{{F{t)F{t') + F{t')F{t))) GS = ±5>fc* ( coth ^) «W*-f) (A10) 

and 

«F(i)» GS = (All) 

Method-III: Our Distribution Function 

In this method we bypass the operator ordering of F(t)s. Instead, we start from coherent state average of the noise 
operator F(t), as given by 10, i.e., the c-number noise 

(F(t)) = Y, [{(<?j(0)) - (i(0))} Kj cosujjt + (p^O))^/ 2 sinwjfl . (A12) 

3 

and then calculate the product of two-time c-numbcrs (F(t))(F(t')). In our notation we have written 

(F(t)){F(t'))=F(t)F(t') (A13) 
Expressing (q'j{0) and in terms of aj and a* as stated earlier, we may write 



F(t)F(t') = |e (^K* 



+ aj) cos ujjt + i\j —^-UJj(a* — aj) sinuijt 



{^ k \j ~^y* k + a k) c °su k t' + i\J ^^uj k (a* k - a fc )sinw fe t'^ j . (A14) 



We now introduce our distribution function (JSJ) for ensemble averaging (...)s over the product F(t)F(t') 



P 3 (a 3 (0),a*(0)) = Nexp 



a, 



n 3 + 1/2 
[ / (v^t 

N exp 



fe-(0)) 2 +^(^.(0)) 2 



2hui 3 (nj + 1/2) 



(A15) 



Noting that 



(af)s = (a* 2 )s = and {otja%) s = {rij + \)8jk, (A16) 



we obtain on statistical averaging over F(t)F(t') 

(F(t)F(t')) s = ( coth ^pf) - f) (A17) 

Several points are now in order. First, we note that the fluctuation-dissipation relation has been expressed in several 
different but equivalent forms in Eas. i|Al|) . ijAlOfl and i|A17(l (or Eq.ljSJ). Because of the prescription of operator 
ordering, the left hand sides of Eas. (|Al(l and (|A10|) do not permit us to identify a c-number noise term. Clearly 
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Ea. (|A17|l suits this purpose and quantum noise as a number F(t) can be identified in (|A17|I or Eq.(|SJ). The ensemble 
averaging procedure of method-Ill is therefore followed in this work. Second, our distribution JSJ is distinctly different 
from the Glauber- Sudarshan distribution l|A8l) because of a vacuum term 1/2 in the width of the distribution (JHJ - 
Third, the origin of distinction of the methods II and III lies in the fact that in the latter case we need not carry 
out an operator ordering of noise as in the former case since we are esssentially dealing with a correlation of numbers 
rather than operators. The uncertainty principle is taken care of through the vacuum term 1/2 in the width of the 
distribution (jSJ) rather than through operator ordering as in method-II. Fourth, the distribution (JSJl as well as the 



Glauber-Sudarshan distribution (|A8|> reduces to classical Maxwell-Boltzmann distribution in the limit fkjj <C fcsT. 
Being a quasi-classical distribution, (|A8|) becomes singular in the vacuum limit n — > 0. The distribution (JSJ on the 
other hand remains a valid non-singular distribution of quantum mean values whose width is simply the natural width 
due to non-thermal vacuum fluctuations. 



APPENDIX B: EVOLUTION EQUATIONS FOR HIGHER-ORDER QUANTUM CORRECTIONS 

The equations for quantum corrections upto fourth order are listed below. 
Equations for the second order are: 



dt 



(5X 2 ) = (5X5P + 5PSX), 

(5P 2 ) = -2T(SP 2 ) - V"(5X5P + 5P5X) - V"'(SXSPSX), 
(5X6P + SPSX) = -T(SXSP + 6P5X)2(6P 2 ) - 2V"(5X 2 ) - V'" (SX 3 ) 



dt 
d 

dt 



Those for the third order are: 



>' s > 



= 3(SXSP8X), 



— (SP 3 ) = -3T(SP 3 ) - W'(SPSXSP) + V" -(SX 2 )(SP 2 ) - -(6PSX 2 5P) + H 2 , 



dt 



dt 
d 

dt 



(SX5PSX) 
(SPSX6P) 



-r(sxspsx) + 2(spsxsp) - v"(sx 3 ) - 

-2T(SP5X5P) + (SP 3 ) - 2V"(5X5P5X) 



V" 



((6X 4 ) - (SX 2 



V' 



((SX 2 ) (SXSP + SPSX) - (SX 3 SP + SPSX 3 )) , 



And lastly, the fourth order equations are: 



(Bl) 



(B2) 



2(SX 3 SP + SPSX 3 ), 

-4r(<5P 4 ) - 2V"(SXSP 3 + SP 3 SX) + 2V"'(SX 2 )(SP 3 ), 

-T(SX 3 SP + SPSX 3 ) - 2V"(SX 4 ) - 'Sh 2 + 6(8PSX 2 SP) 
+ V'" (SX 2 ) (SX 3 ), 

(SXSP 3 + SP 3 SX) = -3T(SXSP 3 + SP 3 SX) + 2(SP 4 ) + W"(h 2 - 2(SP5X 2 SP)) 
+ W'"(SX 2 )(SPSXSP), 

-2T(5PSX 2 5P) - V"(SX 3 SP + SPSX 3 ) + (SP 3 SX + SXSP 3 ) 
V"'(5X 2 )(SXSPSX). 



j t (SX 3 SP + 5PSX 3 ) 



d_ 

dt 



(B3) 



j t (SP8X 2 5P) 
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The derivatives of V(X), i.e., V" or V'" etc. in the above expressions are functions of X. 
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t 

FIG. 1: Plot of correlation function c(t) vs. t (units arbitrary) as given by Eq.(20) (solid line) and Eq.(22) (dotted line) for 
the set of parameter values mentioned in the text. 
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FIG. 2: Plot of inverse of quantum mean first passage time Q M p PT vs. dissipation parameter F (units arbitrary) for different 
temperatures: (a) T = 3.0 (b) T = 1.0 and (c)T = 0.0. 




FIG. 4: Plot of logarithm of inverse of quantum mean first passage time ln( qWFpt ) vs - (units arbitrary) for three 

different V values, V = 2.0, T = 8.0 and T = 20.0. 



